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Thermal  energy  storage  (TES)  has  the  capability  to  absorb,  store  and  release  heat  based  on  dynamic 
surrounding  environmental  conditions.  Sensible  energy  storage  captures  or  releases  heat  with  changes 
in  material's  temperature  while  latent  heat  is  associated  with  a  phase  change  at  an  isotherm  or  near 
isothermal  temperature.  Latent  heat  storage  such  as  using  a  phase  change  material  (PCM)  gains 
growing  attentions  recently  due  to  its  ability  of  storing  significant  thermal  energy  within  a  small 
volume,  making  it  one  of  most  promising  technologies  for  developing  energy  efficient  buildings. 
To  quantify  their  technical  and  economic  feasibility  for  building's  applications,  computational  models 
of  TES  that  can  be  integrated  into  whole  building  energy  simulations  are  highly  demanded.  This  paper 
reviews  the  different  modeling  methods  generally  used  for  PCM  simulations.  A  few  numerical  modeling 
methods  are  observed  in  literature  for  modeling  PCMs  including  the  enthalpy  method,  the  heat  capacity 
method,  the  temperature  transforming  model,  and  the  heat  source  method.  The  study  compares  and 
highlights  the  advantages,  disadvantages  and  limitations  of  these  models  and  methods.  It  particularly 
explores  the  viability  of  these  methods  for  building  applications.  The  paper  further  reviews  the  PCM 
models  that  have  been  integrated  into  prevalent  whole  building  simulation  programs  such  as 
EnergyPlus,  TRNSYS,  ESP-r  etc.  The  study  reveals  that  the  heat  capacity  method  is  mostly  used  in 
programs,  despite  of  its  limitations  on  time  and  spatial  resolutions.  Further  research  is  found  necessary 
to  identify  the  efficiency  and  accuracy  of  these  methods  in  building  applications. 
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1.  Introduction 

Thermal  energy  storage  (TES)  or  thermal  mass  is  a  property  of 
materials  that  describes  its  ability  to  absorb,  store  and  release  heat 
depending  on  the  surrounding  environmental  conditions.  Traditional 
architecture,  for  example,  is  distinguished  with  its  heavy  weight  and 
thermally  massive  construction  elements  to  moderate  the  indoor 
environment  extremes  experienced  in  hot  or  cold  days.  The  thermal 
properties  of  construction  elements  have  significantly  improved 
thermal  comfort  by  manipulating  the  indoor  air  temperature  without 
the  need  of  mechanical  air  conditioning  systems  [1].  On  the  other 
hand,  light  weight  buildings  are  characterised  by  its  lower  thermal 
mass  and  thus  expose  to  significant  temperature  swings,  demanding 
high  cooling  and  heating  energy.  A  dynamic  thermal  mass  such  as 
phase  change  materials  (PCM)  has  been  considered  as  a  promising 
technology  to  reduce  the  inherited  climatic  deficiency  in  light  weight 
buildings. 

The  apparent  advantage  of  using  PCMs  lies  on  the  amount  of 
latent  heat  a  thin  PCM  layer  can  store  compared  to  that  in  a 
sensible  heat  storage  material  such  as  concrete.  For  instance,  a 
wall  of  25  mm  thickness  with  PCM  could  store  an  equivalent 
amount  of  thermal  energy  as  a  420  mm  thick  concrete  wall  [2], 
As  a  result,  the  use  of  PCMs  has  recently  attracted  great  attentions 
for  improving  thermal  and  energy  performance  of  buildings  [3-8]. 
Recent  studies  that  review  PCMs  in  building  applications  can  be 
found  in  [9—1 8].  Various  challenges  are,  however,  arisen  when 
using  PCM  in  buildings,  including  the  variety  of  materials  avail¬ 
able,  material  liability  and  safety,  cost  and  economic  feasibility, 
design  configurations,  integration  with  other  sustainable  energy 
technologies,  impact  on  thermal  and  energy  performance.  The 
problem  can  then  be  considered  as  an  optimization  dilemma 
where  all  these  counterparts’  challenges  are  becoming  critical  in 
the  design  process.  As  a  result,  computational  modeling  is  often 
used  as  an  effective  tool  to  quantitatively  understand  and  help 
resolving  this  optimization  problem. 

This  review  paper  adds  to  the  already  existing  studies  that 
discuss  the  mathematical  modeling  of  phase  change  materials  for 
building  applications  such  as  those  described  in  [19-21].  The 
objective  of  this  study  is  to  systematically  review  the  general 
modeling  theories  and  techniques  for  PCM,  with  an  emphasis  on 
the  specific  models  used  for  simulating  the  thermal  and  energy 
performance  of  PCMs  embedded  in  building  enclosures.  Further¬ 
more,  this  paper  reviews  and  summarizes  the  capabilities, 
limitations  and  validations  of  prevalent  whole  building  simula¬ 
tion  programs  that  have  been  used  for  modeling  phase  change 
materials. 


must  be  met.  For  pure  materials  there  is  a  clear  distinction 
between  the  solid  and  liquid  phase  separated  by  a  sharp  moving 
interface  and  hence  melting  occurs  at  isothermal  temperature. 
For  conduction-dominated  heat  transfer,  the  governing  equation 
can  be  written  for  the  solid  and  liquid  phase,  respectively,  which 
have  to  be  satisfied  by  the  Stefan  condition  as  follows  [22]: 

Heat  transfer  in  the  solid  phase: 
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Heat  transfer  in  the  liquid  phase: 
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The  Stefan  condition  that  enforces  the  heat  balance  at  the 
solid-liquid  interface  is: 
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Very  few  analytical  solutions  are  available  in  the  closed  form 
for  phase  change  problems  and  can  be  found  in  advanced  heat 
transfer  books  such  as  those  by  Crank  [23],  Alexiades  and 
Solomon  [24],  and  Ozis  ik  [25].  Therefore,  approximate  numerical 
solutions  are  usually  used  to  handle  this  class  of  problems.  The 
numerical  methods  for  addressing  these  problems  have  been 
reviewed  in  literature  [26-29]  and  can  be  generally  divided  into: 


1.  Fixed  grid  method  (i.e„  weak  solution):  These  methods  consist 
of  fixed  space  grids  where  the  boundary  is  tracked  by  the  use 
of  an  auxiliary  function.  Different  approaches  are  employed  to 
account  for  latent  heat  evolution  [27,29,30].  This  class  of 
methods  has  been  widely  used  and  therefore  will  be  the  focus 
of  this  paper. 

2.  Deforming  grid  method  or  front  tracking  scheme  (i.e.,  classical 
solution  or  strong  numerical  solution):  These  methods  allow 
the  grid  nodes  move  along  with  the  moving  boundary  layer 
and  thus  the  space  girds  deform  as  the  solution  develops. 
Here  the  interface  is  explicitly  tracked  using  the  Stefan  condi¬ 
tion  [24]. 

3.  Hybrid  method:  These  methods  utilize  the  features  of  both 
fixed  and  deforming  grids  which  uses  a  fixed  background  grid 
and  employs  local  front  tracking  schemes  to  follow  the  move¬ 
ment  of  the  boundary  [26]. 


2.  General  formulation  of  phase  change  problems 


3.  Numerical  formulation  of  phase  change  problems  using 
fixed  grid  methods 


The  main  feature  of  phase  change  problems  (i.e.,  Stefan 
problems)  is  the  moving  boundary  where  the  Stefan  condition 


An  intuitive  approach  in  solving  phase  change  problems  is  to 
explicitly  follow  the  moving  boundary  using  the  front-tracking 
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Nomenclature 

G 

half  range  melting  temperature 

P 

density 

a 

matrix  coefficient 

c 

specific  heat  capacity 

Subscripts  and  superscripts 

f 

fluid  fraction  of  PCM 

h 

enthalpy 

A 

Apparent 

k 

thermal  conductivity 

avg 

average 

L 

latent  heat  of  fusion 

e 

east  node 

n 

the  unit  normal  on  the  phase  interface 

eff 

effective 

s 

source  term 

l 

liquid  state 

T 

temperature 

m 

melting 

t 

time 

n 

iteration  level 

V 

velocity  of  the  interface 

p 

point  node 

X 

space  distance 

s 

solid  state 

w 

west  node 

Greek  letters 

methods.  However,  this  method  needs  to  make  a  priori  assump- 

leads  to  the  following  discretized  equation: 

tion  that  the  boundary  is  smooth  or  monotonic  during  the  period 
[31].  This  assumption  is  not  always  true  and  therefore  reformu- 

K+ 1 

—  /Ip  +  aw+1  x  ^w+1  +ap+1  x  Tp+1  +  a"+1  x  Tne  +  ' 

(5) 

lating  phase  change  problems  using  the  fixed  grid  techniques 
becomes  an  obvious  alternative  [23,30,32,33],  The  Stefan  condi¬ 
tion  Eq.(3)  within  the  fixed  grid  method  is  implicitly  treated  by 
the  reformulated  governing  equation  and  hence  the  position  of 
the  moving  boundary  is  known  when  the  solution  is  converged. 

The  fixed  grid  method  is  simple  compared  to  the  others,  most 
versatile,  convenient,  adaptable  and  easily-programmable  [24], 
The  latent  heat  evolution  is  accounted  for  in  the  governing 
equation  by  using  either  enthalpy  method  [34-38],  heat  capacity 
method  [39-42],  temperature  transforming  model  [43-46],  heat 
source  method  [38,47-50],  or  other  methods  [27,29,51],  The 
following  sections  will  describe  the  widely  used  methods. 

3.3.  The  enthalpy  method 


According  to  Eq.  (5),  it  is  clear  that  the  current  enthalpy  (/ij]+  ) 
is  dependent  on  the  current  value  of  temperature  (TjJ+1)  and 
therefore  the  enthalpy  term  is  nonlinear.  The  equation  cannot  be 
solved  without  using  proper  numerical  techniques  to  handle  this 
nonlinearity.  This  has  to  be  solved  either  by  nonlinear  solvers 
such  as  the  Newton’s  method  or  by  linearizing  the  nonlinear 
terms  and  utilizing  iterative  methods  as  fully  explained  by 
[24,30,31,34,35,52-56].  If  a  non-linear  solver  is  selected,  an 
auxiliary  temperature-enthalpy  function  is  required  for  Eq.  (5) 
and  can  be  written  for  materials  that  change  phases  at  specific 
temperature  range  as  follows  [22]: 


£:>  hp  <  Cs  x  (Tm  —  e ) 

C,x(Tm-e)  <  hp<Csx(Tm+e)+L 


In  the  Enthalpy  method,  the  latent  and  specific  heat  are  com¬ 
bined  into  an  enthalpy  term  in  the  governing  equation.  The  enthalpy 
method  was  proposed  by  Eyres  [38]  to  deal  with  variations  of 
thermal  properties  with  respect  to  temperature.  For  conduction- 
dominated  heat  transfer,  the  governing  Eqs.  (l)-(3)  can  be  reformu¬ 
lated  into  one  equation  where  the  latent  heat  is  absorbed  into  the 
enthalpy  term  as  follows: 


dh  8  f,  8T\ 
P~di  ~  dx\(  8x ) 


(4) 


To  demonstrate  this  method,  a  fully  implicit  control  volume 
approximation  of  Eq.  (4)  for  a  typical  grid  shown  in  Fig.  1 


◄ - ► 

AX 

w 

P 

E 

V 

J 

J 

J 

w 

◄ - - - ► 

i  i 

CD 

\  r 

- - - ► 

5XW  5Xe 


Fig.  1.  A  typical  control  volume  grid. 


hp-(Cs-Ci)xTm-L 

c, 


hp  >  Q  x  (Tm  +  g  )+L 


(6) 


Alexiades  and  Solomon  [24]  have  outlined  numerical  schemes 
for  solving  phase  change  problems  with  the  enthalpy  method 
using  both  linear  and  nonlinear  approaches.  Knoll  on  the  other 
hand  reviewed  various  approaches  utilizing  nonlinear  solvers  to 
resolve  the  Stefan  problem  [55].  He,  in  particular,  developed  an 
algorithm  to  solve  the  Stefan  problem  using  the  Jacobian-free 
Newton-Krylove  method  and  applied  for  two  scenarios:  (1)  pure 
materials  where  melting  occurs  at  isothermal  temperature  and 
(2)  non-isothermal  case  where  phase  change  occurs  at  a  range  of 
melting  temperature. 

An  alternative  approach  to  solving  the  discretized  Eq.  (5)  is 
to  linearize  the  nonlinear  term,  hp+1(T),  using  the  methods 
explained  by  Patankar  [57],  The  discretized  nonlinear  equation 
becomes  linear  with  one  primary  dependent  variable  “Temperature" 
that  can  be  iteratively  solved  with  enthalpy  using  common  linear 
solvers  such  as  direct  methods  (e.g.,  Gauss  elimination  or  tri¬ 
diagonal  algorithm)  or  iterative  methods  (e.g.,  Gauss-Seidel 
method).  Shamsunder  [34],  for  example,  proposed  a  Gauss-Seidel 
iterative  scheme  where  the  solution  sweeps  from  west  to  east  to 
determine  the  state  of  phase  change  and  subsequently  determine 
the  new  nodal  enthalpy.  The  nodal  temperatures  are  then  deter¬ 
mined  based  on  the  discrete  form  of  the  enthalpy-temperature 
relationship.  To  avoid  excessive  iterations,  the  scheme  was  later 
improved  by  introducing  an  over-relaxation  parameter  that  is  used 
at  nodes  where  no  phase  change  occurs  [58].  The  scheme  was 
however  intended  for  phase  change  that  occurs  at  isothermal 
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temperature.  An  iterative  Newton  linearization  scheme  was  intro¬ 
duced  by  Furzeland  [31].  The  solution  process  is  the  same  as  that  of 
Shamsunder  except  that  the  over-relaxation  parameter  can  be 
applied  at  all  nodes. 

Iterative  methods  such  as  Gauss-Seidel  are  inherently  slow 
and  computationally  inefficient.  Therefore,  fast  numerical 
schemes  have  been  introduced  to  improve  the  computational 
efficiency  [35,42,48].  Pham  [42]  proposed  a  method  that  com¬ 
bines  the  features  of  the  enthalpy  and  heat  capacity  methods.  The 
method  consists  of  two  steps:  a  prediction  step  followed  by  a 
correction  step  as  shown  in  Fig.  2.  Based  on  guessed  values,  the 
new  nodal  temperatures  are  predicted  (point  (2)  on  the  graph). 
The  enthalpy  is  determined  based  on  the  predicted  temperature 
values.  The  predicted  temperatures  are  subsequently  corrected  to 
be  consistent  with  the  enthalpy-temperature  curve  (point  (3)  on 
the  graph).  This  temperature  correction  step  is  the  key  of  this 
method.  This  method  is  later  known  to  be  the  "Quasi-Enthalpy" 
method  [59]. 

Voller  pointed  out  that  this  method  might  not  conserve  energy 
at  every  time  step  [22]  and  a  better  conservative  iterative  scheme 
was  proposed  by  Swaminathan  and  Voller  as  illustrated  in  Fig.  3 
[35].  The  method  iterates  the  predicted  and  corrected  intermedi¬ 
ate  values  until  the  convergence  is  achieved.  The  method  has 
been  recently  investigated  as  an  alternative  to  overcome  the 
limitations  of  the  PCM  simulation  algorithm  implemented  in 
ESP-r  [60], 


3.2.  The  heat  capacity  method 

The  heat  capacity  term  in  the  governing  equation  imitates  the 
effect  of  enthalpy  (sensible  and  latent  heat)  by  increasing  the  heat 
capacity  value  during  the  phase  changing  stage.  Two  approaches 
are  generally  used  to  account  for  the  latent  heat  liberation:  the 
apparent  heat  capacity  [22,27,29]  and  the  effective  heat  capacity 
[61,62].  Although  the  two  approaches  differ  in  the  heat  capacity 
approximation,  recent  literatures,  however,  use  the  terminologies 
interchangeably.  More  details  on  the  effective  heat  capacity 
concept  are  explained  by  Poirier  [63]. 

The  apparent  heat  capacity  method  was  introduced  by 
Hashemi  and  Sliepcevich  [64]  to  solve  a  one-dimensional  heat 
transfer  with  phase  change  in  a  mushy  region.  The  conduction- 
dominated  one-dimensional  heat  transfer  equation  using  the 


Fig.  2.  Corrective  non-iterative  scheme  in  the  Quasi-Enthalpy  method  at  a  node 
during  one  time  step. 


Fig.  3.  Corrective  iterative  scheme  in  the  Enthalpy  method  at  a  node  during  one 
time  step. 


apparent  heat  capacity  can  be  written  as: 


p  X  Ca(T)  X 


8T 

8t 


(7) 


The  method  receives  the  popularity  because  the  temperature 
is  the  only  prime  variable  that  needs  to  be  solved  in  the 
discretized  form.  The  key  in  this  approach  lies  in  the  heat  capacity 
approximation.  Two  methods  are  commonly  used  to  approximate 
the  apparent  heat  capacity  term  in  Eq.  (7):  the  analytical/empiri¬ 
cal  relationships  and  the  numerical  approximations. 


3.2.1.  The  analytical/empirical  relationships 

The  heat  capacity  of  a  PCM  can  be  determined  from  the  testing 
data  with  differential  scanning  calorimeters  (DSC).  Manufacturers 
of  PCMs  normally  provide  limited  data  pertained  to  their  products 
such  as  melting  temperature,  heat  of  fusion  and  heat  capacity  at 
solid  and  liquid  states.  Such  minimal  data  can  be  used  to 
approximate  the  heat  capacity  of  a  PCM  using  a  simple  direct 
relationship  with  an  introduction  of  fictitious  melting  tempera¬ 
ture  range  (2  x  e)  [22,43]: 

!CS,  T <Tm- e  (Solid  region) 

Tm-e  <T<Tm+e  (Mushy  region)  (8) 
Q,  T  >Tm  +  e  (Liquid  region) 

Convergence  might  be  an  issue  when  solving  Eq.  (7),  if  the  half 
phase  change  range  (e)  is  set  too  small  or  the  time  step  is  too 
large.  There  is  a  possible  risk  of  missing  the  latent  heat  contribu¬ 
tion  in  a  large  time  step.  Hence,  DSC  testing  results  can  be  used  to 
form  an  empirical  expression  to  approximate  the  heat  capacity. 
Fang  [65],  for  instance,  proposed  a  mathematical  expression  for 
the  heat  capacity  of  paraffin — based  PCM  obtained  from  DSC. 
Others  have  suggested  and  used  alternative  forms  to  approximate 
the  heat  capacity  [66-69]. 


3.2.2.  The  numerical  approximations 

Numerical  approximation  is  an  alternative  when  detailed 
information  about  PCM’s  thermal  behaviors  is  available.  Many 
numerical  approximations  have  been  proposed  in  literature 
[41,70-75].  For  example,  Comini  [70]  applied  a  numerical  tech¬ 
nique  in  the  finite  element  method  where  the  heat  capacity  was 
determined  using  a  derivative  of  enthalpy  with  respect  to  tem¬ 
perature.  Later,  Morgan  [71]  has  improved  the  relationship  to 
avoid  the  convergence  problems.  When  using  an  iterative  scheme, 
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the  heat  capacity  can  be  approximated  using  the  successive 
temperature  and  enthalpy  solutions.  The  temporal  averaging 
proposed  by  Morgan  [71  ]  is  illustrated  in  Fig.  4  and  is  represented 
by  the  following  equation: 


is,  however,  not  commonly  used  but  offers  an  alternative  solution 
when  compared  to  the  apparent  heat  capacity  method. 

3.4.  The  heat  source  method 


A  =  Ah  =  hn-  /t"-1 
AT  Tn-  Tn_1 


(9) 


On  the  other  hand,  Lemmon  [73]  proposed  an  approximation 
based  on  the  space  average  rather  than  the  time  average 
approach.  The  temporal  and  space  average  approximations  are, 
however,  prone  to  convergence  issues  unless  some  precautions 
are  taken  [76].  Solutions  to  the  limitations  of  the  apparent  heat 
capacity  method  have  been  proposed  in  literature  [33,61,77-79]. 
Voller  [22]  found  that  the  apparent  heat  capacity  approximation 
based  on  the  direct  relationships  are  more  accurate  than  the 
Morgan  approximation  used  for  the  cases  he  studied. 


3.3.  The  temperature  transforming  model 


The  temperature  transforming  model  was  developed  by  Cao 
and  Faghri  [80]  to  overcome  the  time  and  spatial  limitations  in 
the  heat  capacity  method.  The  model  has  been  used  by  Faghri  and 
his  co-workers  for  many  applications  [44,45,81].  The  method  is 
also  called  “the  improved  temperature-based  equivalent  heat 
capacity  method”  [82].  While  the  method  was  tested  against 
several  benchmark  examples,  it  has  been  reported  to  produce 
inconsistent  results  especially  when  mass  transfer  through  PCM  is 
considered.  Corrections  were  proposed  to  improve  the  accuracy 
[81,82].  The  key  of  this  method  is  that  the  energy  Eq.  (4)  is 
transformed  into  a  nonlinear  Eq.  (10)  with  a  single  dependent 
variable  “Temperature”  [43]. 
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This  source  term  is  represented  by  the  following  Equation  [43]. 
|Cj  x  e,  T  <  Tm—  e 

S(T)  =  <  (  '  2  Tm—  e<T<Tm+e  (11) 
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The  latent  heat  during  the  phase  change  stage  is  represented 
by  a  source  term  in  the  governing  equation  with  the  heat  capacity 
term  similar  to  the  apparent  heat  capacity  method.  The  method 


Fig.  4.  Apparent  heat  capacity  approximation  at  a  node  during  one  time  step  using 
iterative  methods. 


Using  the  heat  source  method,  the  total  enthalpy  in  the 
governing  Eq.  (4)  is  split  into  the  specific  heat  and  latent  heat 
where  the  latent  heat  acts  as  a  source  term  [23,47].  Eq.  (4)  thus 
becomes: 
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The  method  was  alluded  by  Eyres  [38]  in  the  mid- 1940s. 
In  popular  schemes,  the  phase  change  front  is  tracked  by  the 
evaluation  of  a  nodal  liquid  fraction  field  which  takes  a  value  of 
0  for  solid,  1  for  liquid,  and  a  value  in  the  range  of  [0-1]  for  the 
mushy  region  [23,49],  With  this  approach,  the  fluid  fraction  is 
linearized  and  the  equation  can  be  solved  iteratively  with  tem¬ 
perature.  The  liquid  fraction  can  be  approximated  using  the 
following  auxiliary  equation  [49]: 

fO,  if  T  <  Tm—  e 

fi  =  <  (i,  7',') ,  if  Tm-  e  <  T<  Tm+  e  (13) 

ll,  if  T>  Tm+  e 


When  discretizing  Eq.  (12)  with  a  fully  implicit  scheme  and 
linearizing  the  source  term  “liquid  fraction”  at  the  current  time 
step,  the  discretized  equation  becomes  linear  and  needs  to  be 
solved  for  temperature  in  an  iterative  manner  with  the  liquid 
fraction.  Costa  [50]  has  used  this  method  to  numerically  simulate 
the  latent  heat  thermal  storage. 


3.5.  Summary 

Different  mathematical  models  and  methods  have  been  sug¬ 
gested  in  literature  to  deal  with  phase  change  problems  using  the 
fixed  grids  method:  enthalpy,  heat  capacity,  temperature  trans¬ 
forming  method,  and  heat  source  method.  Every  method  has  its 
main  distinct  feature  for  the  latent  heat  liberation  with  advantages 
and  disadvantages.  Table  1  summarizes  these  methods,  and  high¬ 
lights  the  main  feature  and  their  advantages  and  disadvantages.  For 
many  reasons  including  computational  efficiency,  modeling  accu¬ 
racy  and  flexibility  in  selecting  solution  schemes,  the  enthalpy 
method  is  merited  to  be  an  attractive  mathematical  model  over 
others  for  simulating  phase  change  problems.  In  particular,  it 
becomes  appealing  when  the  corrective  iterative  scheme  (i.e.,  a 
fast  and  energy  conservative  approach),  or  non-iterative 
scheme  (i.e.,  a  quick  but  conservative  approach  at  low  time  steps) 
are  implemented.  To  further  exploit  these  two  features  for  large 
time  steps,  a  quick  but  energy  conservative  approach  is  envisioned. 


4.  Models  for  building  enclosures  with  PCM 

A  few  models  have  been  developed  to  solve  phase  change 
problems  on  the  basis  of  the  general  mathematical  methods 
described  above.  A  list  of  the  models  for  various  engineering  fields 
including  building  applications  was  reviewed  recently  by  these 
studies  [19,20,83-86].  This  section,  however,  provides  a  more  con¬ 
centrated  and  in-depth  review  on  the  models  that  are  proposed  and 
used  for  simulating  PCM  integrated  within  building  enclosures. 

A  few  innovative  and  sustainable  designs  have  been  proposed 
by  integrating  PCM  within  building  construction  elements.  These 
designs  demand  different  levels  of  model  complexity  to  evaluate 
the  thermal  performance  of  such  elements.  The  computational 
models  are  classified  hereafter  as  the  simplified,  intermediate  and 
sophisticated  models.  Within  this  context,  the  simplified  models 
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Table  1 

Feature,  advantages  and  disadvantages  of  mathematical  methods  used  for  phase  change  problems. 


Mathematical 

model  for 

latent  heat 
evolution 

Main  feature 

Advantages 

Disadvantages 

Possible  solution  schemes 

References 

Enthalpy 

Enthalpy  accounts  for 

•  Fast  if  proper 

•  Difficult  to  handle  supercooling 

Iterative  scheme  with  Non-linear 

[24,55] 

method 

sensible  and  latent  heat 

scheme  is  selected 

problems 

solvers  (e.g.,  Newton’s  methods) 

•  Deal  with  sharp  as 

•  The  temperature  at  a  typical  grid 

Linearized-Enthalpy:  Corrective 

[22,35] 

well  as  gradual 

point  may  oscillate  with  time 

iterative  scheme 

phase  change 

Quasi-Enthalpy:  Non-iterative 
Temperature  correction  scheme 

[42] 

Heat  capacity 

Heat  capacity  accounts 

•  Intuitive  since 

•  Lack  of  computational  efficiency 

Iterative  Scheme  (e.g.,  Gauss- 

[22,27,29,61,63,64] 

method 

for  both  sensible  and 

dealing  with  one 

•  Small  time  step  and  fine  grids  are 

Seidel  iterative  scheme  )  if  a 

latent  heat 

dependent 

required  for  accuracy 

proper  heat  capacity 

variable 

•  Difficult  in  handling  cases  where 

approximation  is  selected 

“Temperature” 

the  phase-change  temperature 

•  Easy  to  program 

range  is  small 

•  Suitable  for 

•  Difficult  to  obtain  convergence  with 

gradual 

this  technique,  and  there  is  always 

phase  change 

a  chance  that  the  latent  heat  is 
underestimated 

•  Not  applicable  for  cases  where 

phase  change  occurs  at  fixed 
temperature 

Temperature 

Heat  capacity  and  source 

•  Deal  with  sharp 

•  Not  a  common  method  and 

Iterative  Scheme  (e.g.,  Gauss- 

[43-46,80,81] 

transforming 

term  are  used  to  account 

and  gradual 

therefore  not  tested  to  evaluate  the 

Seidel  iterative  scheme  )  after 

method 

for  sensible  and  latent 

phase  change 

pros  and  cons 

linearizing  the  source  term 

heat 

•  Handle  large  time 

step  and 
course  grids 

Heat  source 

Latent  heat  is  treated  as  a 

•  Intuitive  due  to 

•  Requires  under-relaxation  and 

Iterative  Scheme  (e.g.,  Gauss- 

[23,38,47,49] 

method 

source  term 

separating  the 

therefore  extra  efforts  is  needed  to 

Seidel  iterative  scheme  )  after 

latent  heat  from 

determine  the  optimum 

linearizing  the  source  term 

sensible 

relaxation  factor 

•  Deal  with  sharp 

•  Lack  of  computational  efficiency 

and  gradual 

•  Problems  with  round  off  errors  if 

phase  change 

melting  occurs  over 
temperature  range 

are  rough  approximations  of  the  physics  in  the  phase  change 
process  but  offers  quick  results.  The  intermediate  models  are  a 
tradeoff  between  the  speed  of  the  simplified  models  and  the 
accuracy  and  flexibility  of  the  sophisticated  models.  The  sophis¬ 
ticated  models  are  created  using  well  validated  numerical 
packages  that  offer  a  choice  of  established  and  optimized  numer¬ 
ical  methods.  This  class  offers  a  high  level  of  accuracy  and 
modeling  flexibility  but  is  computationally  expensive. 

4.3.  The  simplified  models 

Detailed  models  for  simulating  PCM  within  building  enclo¬ 
sures  may  capture  more  physics  of  heat  transfer  process.  How¬ 
ever,  simplified  models  are  sometimes  preferred  to  provide  a 
quick  estimation  of  PCM’s  thermal  performance.  Some  simplified 
models  have  been  developed  with  this  intention  [87-91]. 

A  steady-state  analytical  model  for  evaluating  the  benefits  of  PCM 
in  walls  and  roofs  has  been  proposed  by  Kaushik  [87],  The  model 
used  the  heat  capacity  method  to  represent  the  dynamic  thermal 
storage  of  PCM.  The  model  was  utilized  to  analyze  the  dynamic 
thermal  performance  of  a  free  floating  building  with  PCM  embedded 
in  a  south  wall  fapade  [88].  The  result  for  a  typical  mild  winter  day  in 
New-Delhi  showed  that  the  wall  with  PCM  outperformed  that  of  an 
ordinary  wall.  A  rough  model  utilizing  the  heat  capacity  method  was 
developed  to  characterize  the  heat  transfer  process  and  subsequently 
estimate  the  temperature  trend  in  a  PCM  mixed  with  gypsum  plaster 
board  [89].  The  simplified  model  was  able  to  capture  the  overall 


trend  of  air  temperature  in  the  conditioned  room.  Another  simplified 
physical  model  using  the  R-C  network  method  was  developed  and 
validated  for  three  wall  types:  light,  medium  and  heavy  with  shaped- 
stabilized  phase  change  material  [90].  The  model,  however,  had 
to  use  a  genetic  algorithm  to  identify  the  key  model  parameters: 
resistances  and  capacitances  of  the  wall  layers  to  reach  an  optimal 
node  distribution.  When  the  optimal  parameters  are  identified,  the 
model  can  be  used  to  simulate  the  heat  transfer  process  in  a  wall 
unit  that  has  a  PCM  layer.  Although  the  model  is  intended  to  be 
simple,  multiple  procedures  are  necessary  for  practical  applications. 
The  model  was  however  implemented  to  evaluate  the  energy 
performance  of  an  office  building  with  shaped-stabilized  phase 
change  material  embedded  in  a  wall  unit  [91  ]. 

4.2.  The  intermediate  models 

A  variety  of  intermediate  models  using,  respectively,  the 
enthalpy  method,  heat  capacity  method  and  heat  source  methods 
have  been  developed  for  one,  two  and  three  dimensional  cases  for 
building  enclosures. 

4.2.3.  The  enthalpy  method 

In  the  enthalpy  method,  the  enthalpy  may  be  solved  by 
nonlinear  solvers  with  an  auxiliary  function  (e.g.,  temperature- 
enthalpy  relationship)  or  implicitly  in  the  governing  equation 
using  linearization  techniques.  A  theoretical  analysis  based  on  the 
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enthalpy  method  was  presented  in  a  study  evaluating  a  PCM  in 
a  wallboard  for  solar  energy  storage  [92].  A  semi-implicit  Crank- 
Nicolson  method  was  used  for  numerical  discretization,  which 
was  subsequently  solved  using  the  Newton’s  method.  A  more 
sophisticated  two  dimensional  finite  volume  heat  transfer  model 
based  on  the  enthalpy  method  was  developed  and  validated  to 
explore  the  behaviors  of  phase  change  materials  incorporated  into 
building-integrated  photovoltaic  (B1PV)  system  [93,94].  The  heat 
equation  is  solved  with  an  auxiliary  temperature-enthalpy  func¬ 
tion.  The  model  was  utilized  to  perform  an  optimization  study  of 
commercially  available  PCM  products  embedded  into  cavity-wall 
systems  with  different  wall-PCM  configurations  [95],  In  addition 
to  simulating  heat  transfer  process,  the  model  has  the  capability 
to  solve  the  Navier-Stokes  equation  (i.e.,  the  momentum  and 
mass  equations).  The  model  was  expanded  later  to  evaluate  a 
three  dimensional  heat  transfer  process  with  PCM  [96].  It  was 
found  that  the  3D  model  does  not  offer  additional  accuracy  when 
compared  to  the  previously  validated  2D  model.  Another  example 
of  a  validated  finite  element  3D  numerical  model  based  on  the 
enthalpy  method  has  been  suggested  to  simulate  the  PCM  mixed 
with  common  mortars  for  wall  plaster  [97]. 

Using  the  enthalpy  linearization  approach,  a  model  was 
recently  presented  as  an  alternative  method  for  a  PCM  algorithm 
in  ESP-r,  a  whole  building  simulation  program  [60,98].  The 
MATLAB  simulation  environment  was  used  to  develop  a  one 
dimensional  numerical  model  using  a  corrective  iterative  scheme 
proposed  by  Swaminathan  and  Voller  [35]  based  on  the  enthalpy 
linearization.  The  customized  model  in  MATLAB  uses  the  finite 
volume  method  with  a  Crank-Nicholson  scheme  to  produce  a  fair 
comparison  to  ESP-r.  The  model  has  proven  to  be  accurate  and 
fast  when  compared  to  the  ESP-r  results  for  a  BESTEST  Case  600 
model  configured  with  PCM. 


4.2.2.  The  heat  capacity  method 

Phase  change  materials  for  building  applications  such  as 
Paraffin  melt  or  freeze  over  a  temperature  range  compared  to 
pure  materials  where  phase  change  occurs  at  fixed  temperature 
[15-17,85].  This  property  makes  the  heat  capacity  method  an 
attractive  approach  to  simulating  PCM  in  building  applications. 
Utilizing  MATLAB  package,  a  research  group  has  developed  an 
implicit  one  dimensional  finite  difference  model  for  PCM  in  inner 
wallboard,  ceiling  and  floor  with  the  heat  capacity  method  [99]. 
The  discretized  equation  was  solved  using  the  Gauss-Seidel 
iterative  method.  Although  the  lab  experiments  were  limited 
and  simulation  program  was  incomplete  at  that  stage,  the  overall 
benefits  from  PCM  in  wallboard  were  evident. 

A  semi-implicit  one  dimensional  finite  volume  heat  transfer 
model  for  simulating  PCM  in  a  ceiling  of  a  room  using  the  heat 
capacity  method  was  developed  and  validated  by  Pasupathy 
[100,101].  The  model  was  solved  using  the  tri-diagonal  matrix 
algorithm  (TDMA)  with  very  small  time  step.  Although  the  overall 
trend  of  indoor  air  temperature  was  captured  by  the  model,  the 
numerical  results  were  not  in  a  good  agreement  with  the  experi¬ 
ments  due  to  many  limitations  of  the  model.  The  same  model  was 
later  used  for  evaluating  the  PCM  integrated  into  a  roof  system  [102]. 

The  heat  capacity  method  has  also  been  implemented  in  a  one 
dimensional  numerical  model  to  evaluate  shape-stabilized  phase 
change  materials  embedded  with  a  floor  heating  system  [103]. 
The  specific  heat  capacity  was  used  to  account  for  the  enthalpy 
of  PCM  at  different  temperature  regimes.  The  model  gave  good 
agreement  results  when  compared  to  experimental  data.  The 
model  was  also  used  for  PCM  evaluations  under  different  climates 
and  various  system  configurations  [5,104,105],  A  variety  of 
modeling  applications  of  PCM  embedded  in  floor  system  for 


different  purposes  using  the  heat  capacity  method  have  been 
reported  in  literature  [106-108]. 

PCMs  have  also  been  integrated  in  transparent  building 
envelopes  such  as  glazed  windows.  An  explicit  one  dimensional 
finite-difference  model  based  on  the  heat  capacity  method  was 
extended  to  evaluate  PCM  performance  when  integrated  into  a 
double  glazing  system  [109,110].  The  developed  model  was 
validated  with  experimental  data  and  then  subsequently  utilized 
to  evaluate  the  impact  of  PCM  on  heat  loss  and  gains. 

While  models  developed  for  building  envelope  are  normally  for 
one  dimensional  geometry,  two  and  three  dimensional  heat 
transfer  approaches  have  also  been  suggested  for  modeling  PCM 
using  the  heat  capacity  method.  In  the  early  1990s,  a  numerical 
code  “WALL88”  was  proposed  for  modeling  two  dimensional 
transient  thermal  transport  and  storage  of  both  sensible  and  latent 
heat  [111].  The  model  was  validated  against  analytical  solution 
and  experimental  lab  results.  The  model  was  found  to  give  an 
excellent  agreement  with  experimental  results  only  after  allowing 
the  PCM  to  melt  over  a  temperature  range  rather  than  at 
isothermal  temperature.  A  three  dimensional  finite-difference 
heat  transfer  model  using  the  heat  capacity  method  was  devel¬ 
oped  to  study  the  thermal  performance  of  randomly  mixed  PCM 
and  laminated  PCM-wallboard  systems  [112,113].  Although  the 
numerical  model  was  not  validated  in  these  papers,  the  simulation 
results  were  helpful  to  conclude  that  laminated  PCM-wallboard 
performs  thermally  better  than  the  randomly  PCM-wallboard.  This 
model  was  later  validated  against  an  experiment  and  found  a 
maximum  of  3%  deviation  from  the  average  experimental  results 
[114].  The  model  was  further  used  to  evaluate  the  PCM  applica¬ 
tions  in  the  drywall  in  a  passive  solar  building  [115].  The  results 
confirmed  the  conclusions  from  previous  studies  for  the  applica¬ 
tion  of  laminated  PCM-wallboard.  Optimizing  the  PCM  distribu¬ 
tion  within  building  envelope  is  the  overall  goal  of  an  energy 
efficient  design.  The  heat  capacity  method  was  adopted  by  an  in- 
house  software  “CODYMUR”,  developed  by  a  team  from  France,  to 
optimize  the  use  of  a  PCM  wallboard  for  building  energy  use  [116]. 


4.2.3.  The  heat  source  method 

The  heat  source  method  is  an  intuitive  approach  due  to  the 
separation  of  specific  and  latent  heat.  An  explicit  one  dimensional 
finite-difference  heat  transfer  model  for  a  wall  with  PCM  was 
developed  using  the  heat  source  method  by  Athienitis  [117]  and 
validated  against  a  full-scale  outdoor  test-room  with  PCM  gyp¬ 
sum  board  at  interior  side.  The  model  showed  a  reasonable 
agreement  with  experimental  results.  A  heat  transfer  model  of  a 
newly  developed  hybrid  thermal  energy  storage  system  (HTESS) 
using  PCM  capsules  in  a  wall-unit  was  developed  and  validated 
for  managing  solar  and  electric  energy  [118].  The  numerical 
model  uses  the  heat  source  method  similar  to  that  proposed  by 
Voller  [48]  where  the  latent  heat  evolution  is  represented  by  a 
source  term  in  the  governing  equation.  The  fluid  fraction  is  the 
key  to  track  the  latent  heat  process. 

Phase  change  materials  incorporated  within  floor  systems  was 
evaluated  using  a  one  dimensional  finite  volume  heat  transfer 
model  based  on  the  heat  source  method  [119].  The  numerical 
model  was  validated  with  a  benchmark  analytical  solution  of 
Stefan  problem  explained  by  Hu  [29].  After  optimizing  the  grid 
and  the  time  resolution,  the  developed  model  together  with  an 
optimization  algorithm  was  used  to  perform  an  optimization 
analysis  on  PCM-floor  designs.  A  two  dimensional  numerical 
model  was  developed  based  on  the  heat  source  method  to 
simulate  the  effect  of  PCM  in  the  design  of  a  solar  passive  wall 
[120].  The  model  has  been  verified  using  benchmark  cases 
documented  in  literature.  The  model  was  later  validated  using 
experimental  data  performed  in  the  Lab  and  found  to  be 
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unsatisfactory  due  to  the  limitation  of  handling  the  super-cooling 
effect  inherited  in  the  tested  PCM  [121]. 

4.3.  The  sophisticated  models 

Developing  a  numerical  model  in  two  or  three  dimensional  domain 
is  complicated  and  difficult  to  be  generalized  for  different  geometries, 
applications  and  physics;  hence  existing  simulation  packages  such  as 
COMSOL  (formerly  known  as  FEMLAB)  [122],  ANSYS-FLUENT  [123], 
HEATING  [124]  and  others  are  used  as  convenient  design  tools. 
Although,  these  models  offer  high  level  of  flexibility,  they  are  not 
fully  explored  for  heat  transfer  process  with  phase  change. 

One  study  has  used  a  commercial  package  FEMLAB  (later 
COMSOL  multiphysics)  to  develop  a  wall  model  with  phase 
change  materials  using  the  enthalpy  method  and  heat  capacity 
method  [125],  COMSOL  is  a  finite  element  simulation  package 
that  allows  multi-physics  modeling  for  many  engineering  appli¬ 
cations.  Utilizing  this  package,  the  created  model  was  validated 
against  experimental  results.  It  was  found  that  both  numerical 
methods  give  good  estimation  of  latent  heat  evolution  process. 
However,  the  heat  capacity  method  found  to  be  more  precise  with 
the  experimental  results  when  a  narrow  melting  temperature 
range  of  2  °C  was  selected.  COMSOL  has  also  been  used  to  study 
envelope  systems  with  PCM  [126].  The  numerical  results  from 
COMSOL  were  successfully  compared  with  another  well- 
established  numerical  model  “WUFI-5”.  COMSOL  is  flexible  in 
modeling  multi-physics  within  irregular  and  complex  geometries. 
For  example,  an  innovative  honeycomb  wallboards  with  PCM 
have  been  modeled  in  a  3D  domain  using  the  heat  capacity 
method  [127].  The  simulation  results  showed  a  very  good  agree¬ 
ment  with  the  experimental  results. 


A  computational  fluid  dynamics  (CFD)  simulation  package 
“FLUENT”  has  been  utilized  to  evaluate  a  heat  source  method 
when  using  user-defined  functions  for  heating  and  cooling  cycles 
of  PCM  rather  than  using  one  idealized  function  to  represent  both 
phenomena  [128],  The  results  of  a  PCM  box  model  utilizing  the 
two  functions  showed  very  small  error  (in  root  mean  square 
(RMS)  values)  when  compared  with  a  case  of  using  an  ideal 
function  for  phase  change.  Another  example  of  using  FLUENT  is 
reported  for  PCM  integration  into  a  wall  cavity  system  [129]. 

Heat  Engineering  And  Transfer  In  Nine  Geometries  (HEATING) 
is  a  multidimensional,  general-purpose  heat  transfer  code  that 
has  been  extensively  validated  under  ASHRAE  project  RP-1145 
[130],  The  code  can  also  be  used  to  model  the  phase  change  using 
the  heat  capacity  method.  Ahmad  [131]  has  used  the  program  to 
study  the  behaviors  of  PCM  in  wallboard  of  a  test  cell.  The  model 
was  validated  using  experimental  test  results  and  found  to  agree 
well  with  experiments.  The  PCM  research  program  at  Oak  Ridge 
National  Laboratory  (ORNL)  has  used  this  program  to  study  the 
thermal  behaviors  of  PCM  in  complex  two  and  three  geometries 
in  building  envelope  [132].  Lab  tests  using  a  heat  flow  meter 
apparatus  (HFMA)  have  been  conducted  to  validate  the  model  in 
HEATING.  HEATING  has  also  been  used  as  a  standard  benchmark 
numerical  package  to  validate  the  finite-difference  algorithm 
used  for  PCM  modeling  in  EnergyPlus  [133]. 

4.4.  Summary 

A  variety  of  models  for  different  building  enclosures  have  been 
developed  using  various  simple,  intermediate  and  sophisticated 
approaches.  Table  2  summarizes  the  models,  application  usages, 
and  validations.  It  is  obvious  that  very  few  simplified  models  have 


Table  2 

Modeling  approaches  for  latent  heat  evolution  in  building  enclosure. 


Complexity  Latent  heat  evolution’s  approach 

Building  enclosure 

Modeling 

Solution  strategy 

Validation 

References 

level 

case  studied 

formulation 

Simplified  models 

Heat  capacity  method 

Wall  and  roof 

Steady-state 
analytical  model 

[87,88] 

Wallboard 

Experimental 

[89] 

Optimum  nodes  for  heat  capacity 
distribution  using  genetic  algorithm 

Wall 

R-C  Network 

N/A 

[90,91] 

Intermediate  models 

Enthalpy  method 

Wall 

FVM:  ID 

Newton’s  method 

[92] 

BIPV 

FVM:  2D  &  3D 

Non-linear  solver 

Experimental 

[93,95,96] 

Wall 

FVM:  ID 

Iterative  corrective 
scheme 

Comparative 

[60,98] 

Heat  capacity  method 

Wall 

FDM:  ID 

G-S 

Experimental 

[99] 

Ceiling/roof 

FVM:  ID 

TDMA 

Experimental 

[100-102] 

Floor 

FDM:  ID 

G-S 

Experimental 

[5,103-105] 

Glazed-Windows 

FDM:  ID 

EM 

Experimental 

[109,110] 

Wallboard 

FDM:2D 

Analytical  and 
experimental 

[111] 

Wallboard 

FDM:  3D 

Experimental 

[112-115] 

Heat  source  method 

Wall 

FDM:  ID 

EM 

Experimental 

[117] 

Wall 

FDM:  ID 

TDMA 

Experimental 

[118] 

Wall 

FVM:  ID 

Iterative  scheme 

Analytical 

[119] 

Wall 

FVM:  2D 

TDMA 

Analytical,  comparative 
and  experimental 

[120] 

Sophisticated  numerical  packages 

COMSOL  Heat  capacity  method  and  heat  source 

Wall 

FEM:  2D 

Experimental 

[125] 

method 

Heat  capacity  method 

Wall 

FEM 

Comparative 

[126] 

Heat  capacity  method 

Wall 

FEM:  3D 

Experimental 

[127] 

FLUENT  Heat  source  method 

Wall 

FVM:  3D 

SIMPLE  algorithm 

Experimental 

[128] 

Heat  source  method 

Wall 

FVM:  2D 

Experimental 

[129] 

HEATING  Heat  capacity  method 

Wall  and  roof 

FDM:  ID,  2D  and 

Point-successive  over- 

Experimental 

[131-133] 

3D 

relaxation  iteration 

FVM:  finite  volume  method,  FDM:  finite  difference  method,  FEM:  finite  element  method,  G-S:  Gauss-Seidel  iterative  method,  TDMA:  tridiagonal  matrix  algorithm, 
EM:  explicit  time  stepping  marching,  R-C:  resistance-capacitance. 
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been  suggested  due  to  the  complexity  of  approximating  the  heat 
transfer  process  associated  with  phase  change.  The  intermediate 
models  are  commonly  used  but  are  developed  for  specific  appli¬ 
cations  and  to  investigate  explicit  envelope  designs.  Hence,  they 
lack  flexibility  in  analyzing  complex  and  advanced  design  alter¬ 
natives  which  becomes  norms  for  selecting  optimal  or  near 
optimal  designs.  Sophisticated  models  offer  flexibility  in  solving 
complex  and  multi-physics  problems  but  are  not  fully  explored 
for  modeling  PCMs.  This  is  partly  due  to  the  computational 
inefficiency.  They  additionally  demand  considerable  amount  of 
detailed  data  inputs,  lengthy  model  setup  and  validations,  and 
limited  access  to  the  source  codes. 

Generally,  all  models  that  adopt  the  heat  capacity  or  heat 
source  methods  must  be,  however,  used  with  small  time  steps  to 
attain  acceptable  accuracy  and  therefore  slow  for  whole  year 
simulation  which  is  typical  for  building’s  thermal  performance 
evaluation.  In  addition,  many  existing  models  ignore  inherited 
characteristics  of  some  PCMs  such  as  hysteresis  or  subcooling  and 
therefore  cannot  be  used  for  this  particular  application. 

5.  PCM  models  integrated  into  whole  building  simulation 
programs 

Many  detailed  simulation  programs  are  nowadays  available  to 
assist  designers,  researchers,  manufacturing  companies  to  imple¬ 
ment  new  technologies  and  evaluate  innovative  ideas  that 
improve  the  energy  and  thermal  performance  of  buildings. 
Detailed  simulation  tools  perform  computations  on  an  hourly  or 
sub-hourly  bases  for  accurate  considerations  of  the  dynamic 
interactions  between  all  thermal-based  elements  associated  with 
comfort  and  energy  consumption,  including  building  envelope, 
HVAC  systems,  lighting  and  control  devices  [134].  Many  building 
simulation  tools  are  listed  at  the  U.S.  Department  of  Energy  (DOE) 
web  directory  [135].  The  twenty  prevalent  whole  building  energy 
simulation  programs  that  are  considered  accurate  and  capable  of 
handling  the  dynamic  behaviors  of  a  building  and  its  systems  are 
reviewed  by  Crawley  et  al.  [136].  Few  whole  building  simulation 
programs  can  handle  the  thermal  performance  of  building  envel¬ 
ope  with  phase  change  materials  such  as  EnergyPlus,  TRNSYS, 
ESP-r,  and  BSim.  In  addition,  some  other  programs  with  limited 
capabilities  are  available  for  modeling  phase  change  in  buildings. 
The  following  paragraphs  brief  and  compare  the  conditions  of 
these  programs. 

5.3.  EnergyPlus 

EnergyPlus  uses  the  Conduction  Transfer  Functions  (CTF)  to 
approximate  heat  transfer  in  building  envelope.  Since  the  CTF 
method  uses  the  historical  values  of  heat  flux  in  the  computation, 
Barbour  [137]  has  studied  the  possibility  of  using  this  method  in 
EnergyPlus  to  approximate  the  latent  heat  evolution  in  building 
envelope.  The  study  developed  multiple  sets  of  CTFs  based  on  the 
temperature  of  phase  change  materials.  A  switching  mechanism 
was  proposed  to  exchange  between  these  sets  during  the  simula¬ 
tion.  The  CTF-switching  algorithm  was  found  to  be  within  20% 
accuracy  for  a  range  of  conditions  typically  encountered  in 
buildings. 

The  capability  of  modeling  PCMs  has  been  facilitated  in 
EnergyPlus  program  Version  2.0  released  in  April  2007  by  adding 
a  conduction  finite  difference  (CondFD)  solution  algorithm  [138]. 
The  algorithm  uses  a  semi-implicit  finite  difference  scheme  based 
on  the  heat  capacity  method  with  an  auxiliary  enthalpy- 
temperature  dataset  to  account  for  latent  heat  evolution  [139]. 
Using  this  dataset,  the  heat  capacity  is  approximated  using  a 
temporal  averaging  approach  similar  to  that  proposed  by  Morgan 


[71],  While  the  previous  versions  of  EnergyPlus  had  a  semi- 
implicit  scheme  for  modeling  PCMs,  a  fully  implicit  scheme  has 
been  recently  added  to  Version  7  of  the  program  with  more 
numerical  flexibility  [140].  For  both  schemes,  it  is  however 
recommended  to  use  a  small  time  step  for  accurate  results. 

Experimental  validations  have  been  conducted  for  this  algo¬ 
rithm  with  mixed  feelings  of  accuracy.  Castell,  for  example  [141], 
found  that  EnergyPlus  simulation  results  did  not  show  a  good 
agreement  with  the  experiments  when  phase  change  materials 
were  implemented  in  concrete  blocks.  The  study  concluded  that 
the  simulation  results  did  not  reflect  the  thermal  improvement  of 
PCM  observed  in  the  test  cells.  However,  the  study  highlighted 
that  the  weather  data  used  in  this  simulation  was  not  represen¬ 
tative  of  the  actual  weather  data. 

An  experimental  test  shed  with  a  commercial  PCM  product  has 
been  used  to  validate  EnergyPlus  (Version  5)  simulation  results 
under  the  climatic  conditions  of  Phoenix,  Arizona  [142],  It  was 
found  that  the  predicted  energy  consumptions  were  half  during 
winter  and  slightly  greater  for  the  summer  months.  In  addition, 
the  time  shift  was  observed  for  a  very  short  time  span  during  the 
month  of  April  (3  min)  and  October  (9  min). 

Under  the  climatic  conditions  of  Auckland  New-Zealand,  an 
experimental  study  using  PCM  in  gypsum  board  has  been  com¬ 
pared  to  the  simulation  results  of  EnergyPlus  using  both  historical 
weather  data  and  actual  measured  weather  data  [143].  Although 
EnergyPlus  model  using  actual  weather  data  has  captured  the 
overall  trend  of  indoor  air  temperature  but  failed  to  accurately 
predict  the  actual  indoor  air  temperature  from  measurements. 
The  study  highlighted  that  due  to  many  parameters  including  air 
infiltration,  the  simulation  results  might  deviate  from  the  actual 
measurement. 

An  early  successful  validation  of  the  CondFD  solution  algo¬ 
rithm  used  for  PCM  modeling  has  been  reported  by  Zhuang  [144] 
using  two  envelope  systems  with  PCM:  envelope  “A”  (one  layer  of 
PCM:  melting  temperature  at  40  °C)  and  envelope  “B”  (two  layers 
of  PCM:  one  melting  temperature  at  40  C  and  another  at  33  °C). 
The  study  shows  that  the  largest  relative  difference  in  indoor 
temperature  is  12.41%  and  the  least  is  0.71%  between  the 
simulation  and  testing  results  in  a  sequential  36  h  period  on 
envelope  “A”  condition.  For  envelope  "B"  condition,  the  largest 
relative  difference  is  8.33%  and  the  least  is  0.33%  in  a  sequential 
72  h.  It  was  concluded  that  the  most  important  factors  in  redu¬ 
cing  the  discrepancies  between  the  simulation  and  the  test  results 
are  to  use  proper  actual  weather  data  as  well  as  using  proper 
material  thermal  characteristics.  Other  successful  validations  of 
EnergyPlus  algorithm  for  PCM  have  been  conducted  by  Campbell 
[145]  and  Chan  [146]  using  published  experimental  data  by 
Kuznik  [147],  For  both  validation  studies,  the  indoor  air  tempera¬ 
ture  was  found  to  agree  well  with  the  experimental  results. 

A  field  test  house  was  used  to  study  the  impacts  of  PCMs  in 
building  envelope  and  consequently  used  to  validate  the  CondFD 
algorithm  in  EnergyPlus  [148].  The  test  house  used  the  cellulose 
insulation  mixed  with  20%  PCM  by  mass.  The  study  reported  that 
simulated  daily  average  heat  flux  through  walls  was  within  9%  of 
the  field  measurements.  In  addition,  simulation  results  for  tem¬ 
perature  distribution  through  envelope  compared  fairly  well  with 
the  experimental  data  apart  from  some  delayed  response  com¬ 
pared  to  the  measurement.  However,  EnergyPlus  has  given 
unreasonable  results  for  heat  fluxes  and  temperature  distribu¬ 
tions  in  the  attic  floor  of  the  experimental  house. 

In  addition  to  the  validations  above,  the  EnergyPlus’s  devel¬ 
oper  team  has  performed  rigorous  validation  and  verification 
studies  for  general  heat  transfer  calculations  as  well  as  the 
CondFD  solution  algorithm  [133,149,150],  These  validation  stu¬ 
dies  used  analytical  benchmark  solutions,  comparative  tests  with 
well-established  program  “HEATING",  and  experimental  results. 
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The  studies  concluded  that  versions  prior  Version  7  contains  two 
bugs  and  subsequently  will  be  fixed  in  a  later  version.  It  is  also 
recommended  to  follow  guidelines  and  bear  in  mind  limitations 
in  using  Version  7: 


•  The  time  step  should  be  shorter  than  3  min. 

•  For  an  accurate  hourly  thermal  performance,  1/3  of  the  default 
node  space  should  be  used. 

•  Hysteresis  in  PCM  is  not  modeled  in  EnergyPlus  and  therefore 
inaccurate  results  may  be  produced. 


5.2.  TRNSYS 

TRNSYS  is  a  modular  program  where  components  modules 
“TYPES”  are  linked  together  in  which  output  of  one  type  can  be 
an  input  to  another  in  the  model.  It  has  been  widely  used  for 
modeling  building  and  its  complex  systems.  Due  to  its  modular¬ 
ity,  users  can  either  utilize  available  types  in  the  simulation 
package  or  develop  new  modules  and  easily  integrate  to  the 
TRNSYS  simulation  package.  Many  features  have  been  introduced 
in  Version  16,  including  a  graphical  user  interface  “Simulation 
studio”  and  the  possibility  to  call  external  programs  such  as 
MATLAB,  FLUENT  and  many  others  [151], 

Many  models  have  been  proposed  in  TRNSYS  for  modeling 
phase  change  heat  transfer  in  building  envelope  but  majority  are 
proprietary  research  modules.  Ghoneim  for  example  used  a 
modified  type  of  the  thermal  storage  wall  “TYPE36”  where  the 
use  of  PCM  has  been  tested  for  thermal  storage  in  a  wall  system 
[152,153],  The  model  was  based  on  the  enthalpy  method  and 
solved  using  an  explicit  scheme.  Despite  the  numerical  problems 
encountered  when  modeling  PCM  due  to  the  smaller  time  step 
required  for  the  stability  of  the  explicit  scheme,  the  model  was 
successfully  integrated  into  TRNSYS  and  validated  against  pub¬ 
lished  data  for  a  concrete  storage  wall.  Another  explicit  numerical 
scheme  using  the  enthalpy  method  was  developed  for  modeling 
the  effects  of  integrating  PCM  into  a  solar  wall  [154,155],  The 
module  “TYPE58”  was  integrated  into  TRNSYS  to  explore 
the  significance  of  heating  outside  air  for  ventilation  purposes 
in  the  experimental  house. 

Modeling  PCM  in  TRNSYS  is  recently  provided  through 
“TYPE204"  by  a  team  of  researchers  from  Helsinki  University  of 
Technology  [156],  The  model  simulates  heat  transfer  through  a 
wall  in  a  three  dimensional  domain  using  the  Crank-Nicolson 
scheme  with  729  nodes.  The  model  can  indeed  use  a  fully  implicit 
or  explicit  scheme  with  an  appropriate  selection  of  a  parameter 
that  switches  between  different  schemes.  The  model  uses  the 
heat  capacity  to  account  for  latent  heat  evolution  in  the  wall. 
Although  this  type  has  not  been  validated  in  its  3D  form  due  to  its 
poor  computational  efficiency,  Ahmad  [157]  has  converted  the  3D 
model  into  a  ID  module  “TYPE101”  and  validated  the  modified 
code.  The  simulation  results  were  compared  to  experimental 
results  from  two  test  cells:  one  without  PCM  and  another  with 
PCM.  While  the  model  without  PCM  works  well  when  compared 
to  the  experimental  results,  the  model  with  PCM  overestimates 
the  daily  peak  indoor  temperature  in  the  cell.  The  authors 
outlined  several  reasons  for  this  discrepancy  including:  (i)  eva¬ 
luation  of  the  energy  transmitted  through  the  window,  (ii) 
imprecision  in  the  melting  temperature  range  taken  in  the  heat 
capacity  definition,  (iii)  values  of  the  convective  heat  transfer 
coefficient  between  wall  surfaces  and  ambient  air  and  (iv) 
existence  of  cold  bridges.  Out  of  these  reasons,  it  was  found  that 
correcting  cold  bridges  by  introducing  extra  term  for  resistance 
improved  the  simulation  results  significantly  when  compared  to 
the  experimental  results. 


A  study  reported  a  simplified  approach  of  simulating  PCM  in 
walls/ceiling  and  floor  in  TRNSYS  [158].  The  approach  is  to  use  the 
existing  capability  of  TRSNYS  to  simulate  a  standard  active  wall  in 
“TYPE56”  (i.e„  building  module  in  TRNSYS).  The  key  in  this 
approach  is  a  user  input  of  equivalent  heat  transfer  coefficients 
introduced  in  each  time  step  of  the  simulation  that  characterizes 
the  thermal  behaviors  of  the  wall  with  PCM.  The  model  does  not 
evaluate  the  real  heat  transfer  behaviors  in  PCM  but  accurate 
enough  for  modeling  PCM  thermal  behaviors.  The  model  has  been 
validated  under  laboratory  setting  conditions. 

On  the  other  hand,  Schranzhofer  et  al.  [159]  have  developed  a 
PCM  module  “TYPE241  ”  where  PCM  was  modeled  as  an  internal 
layer  based  on  the  heat  source  method.  TRNSYS  capability  was 
utilized  to  model  other  envelope  layers  using  the  transfer  func¬ 
tion  method  by  creating  dummy  contact  zones  between  the  PCM 
layer  and  the  remaining  layers.  In  this  type,  the  PCM  is  modeled 
using  external  code  based  on  finite  different  method  with  other 
layers  modeled  through  CTF  algorithm  available  internally  in 
TRNSYS  “TYPE56”.  One  advantage  of  this  approach  is  the  short 
computational  time  needed  for  numerical  solution  but  the  phy¬ 
sics  might  not  be  captured  well  because  of  assumptions  involved 
in  the  dummy  contact  zones.  The  model  however  was  not 
validated  due  to  a  lack  of  appropriate  experimental  data. 

Kuznik  et  al.  have  recently  developed  a  new  model  “TYPE260’  in 
TRNSYS  utilizing  the  heat  capacity  method  [160].  The  model  is  semi- 
implicit  since  the  physical  properties  of  PCM  used  in  the  computa¬ 
tions  are  calculated  from  previous  time  step.  This  type  has  been 
validated  with  two  lab  tests  conducted  by  authors:  one  when  the 
outdoor  temperature  was  increased  in  two  steps  and  the  second 
when  it  was  a  sinusoidal  behavior.  The  heating  heat  capacity  curve 
was  used  for  the  numerical  modeling.  For  both  validations,  the 
simulation  results  showed  good  agreement  with  the  test  results. 

A  newly  developed  one-dimensional  heat  transfer  model  using 
the  heat  capacity  method  was  applied  to  a  dividing  wall  with  16 
glass  bricks  filled  with  PCM  in  TRNSYS  [161].  The  model  has  been 
validated  and  showed  fair  agreement  with  experimental  results.  In 
addition,  a  simplified  PCM  module  “TYPE1270’’  has  been  recently 
developed  by  Thermal  Energy  System  Specialists  (TESS)  and  added 
to  its  commercially  available  individual  components  [162].  The 
module  simulates  PCM  as  an  internal  layer  within  an  envelope 
system.  The  model  is  currently  limited  to  materials  that  melt/freeze 
at  isothermal  temperature  and  with  constant  specific  heat  at  solid 
and  liquid.  In  the  transition  state,  the  PCM  layer  temperature  is 
constant  and  the  model  tracks  the  energy  absorption  and  release. 
The  tracking  methodology  is  similar  to  the  heat  source  method  and 
therefore  can  be  identified  as  “Quasi-Heat  Source  Method". 

5.3.  ESP-r 

ESP-r  is  a  dynamic  energy  simulation  tool  of  UK,  used  for 
modeling  thermal,  visual  and  acoustic  performance  of  buildings 
[163].  With  many  features  suitable  to  model  advance  sustainable 
energy  technologies,  ESP-r  has  the  capability  to  model  phase 
change  materials  using  two  methods:  the  effective  heat  capacity 
method  and  the  additional  heat  source  method  [62,164,165],  ESP- 
r  uses  four  models  for  PCM  simulation,  with  one  that  accounts  for 
sub-cooling,  using  special  materials  function.  However,  it  is 
necessary  to  use  a  small  time  step  to  obtain  accurate  results  for 
these  two  methods.  While  simulation  results  using  ESP-r  have 
been  found  in  literature,  none  showed  any  substantial  validations 
for  these  two  algorithms  in  ESP-r  [166-169]. 

5.4.  BSim 

BSim  is  a  dynamic  simulation  program  originated  from 
Denmark  that  offers  an  easy  user  graphical  interface  [170].  Using 
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the  quasi-steady  state  in  building  modeling,  the  program  models 
phase  change  using  the  heat  capacity  method  [171].  The  simula¬ 
tion  time  step  has  to  be  small,  too,  for  accurate  prediction.  Lab 
test  results  from  literature  were  used  to  validate  the  model  on 
three  cases:  continuous  heating,  continuous  cooling,  and  heating 
but  with  initial  temperature  below  melting  point  of  PCM.  The 
simulation  model  captures  the  overall  trend  of  actual  thermal 
behaviors  of  PCM  but  with  small  deviations. 


5.5.  Other  building  simulation  programs 

Some  other  building  simulation  programs  have  been  devel¬ 
oped  for  specific  research  purposes  of  modeling  phase  change 
materials,  such  as  RADCOOL  [172],  ESim  [173],  and  CoDyBa  [66]. 
Few  have  been  proposed  and  developed  to  simulate  simple 
building  configurations,  such  as  PCMExpress  [174]  or  the  one 
using  Engineering  Equation  Solver  (EES)  [106].  Due  to  limited 


Table  3 

Numerical  methods  for  latent  heat  evolution  in  building  simulation  programs. 


Building 

Module 

Numerical 

Numerical 

Time  stepping 

Limitations/Constrains 

Validation 

Reference 

simulation 

identification 

formulation 

method  used  for 

latent  heat 
evolution 

scheme 

EnergyPlus 

CondFD 

FDM:  ID 

Heat  capacity 

1.  Implicit 

•  Time  step  <  3  min 

Analytical, 

[133,144-146,148-150) 

method 

2.  Semi-implicit 

•  Small  grids 

Comparative  and 

•  Hysteresis  in  PCM  is  not 

Experimental 

modeled 

•  Phase  Change  at  isothermal 

temperature  is  not  modeled 

TRNSYS 

Modified 

FDM:  ID 

Enthalpy  method 

Explicit 

•  Low  time  step 

Limited  validation 

[152,153] 

“TYPE36” 

•  No  access  to  the  code 

using 

experimental 
results  for 

concrete 

“TYPE58” 

FDM:  2D 

Enthalpy  method 

Explicit 

No  access  to  the  code 

Experimental 

[154] 

“TYPE204” 

FDM:  3D 

Heat  capacity 

Select  an  appropriate 

Computationally  inefficient 

N/A 

[156] 

method 

factor  for  implicit, 
semi-implicit  or 
explicit 

“TYPE101” 

FDM:  ID 

Heat  Capacity 

Semi-implicit 

•  A  correction  factor  to 

Experimental 

[157] 

Method 

(Crank-Nicolson) 

account  for  cold  bridges  has 
to  be  used  for  model 

accuracy 

TRNSYS 

Equivalent 

Variable  heat 

•  Real  heat  transfer  physics  in 

Experimental 

[158] 

“Active  Wall” 

heat  transfer 

source  function 

PCM  is  not  modeled 

coefficients 

mimicking  PCM 
behavior 

“TYPE241  ” 

FDM:  ID 

Heat  source 
method 

No  Published  data 

N/A 

[159] 

“TYPE260” 

FDM:  ID 

Heat  capacity 

Implicit 

•  Thermal  properties  including 

Experimental 

[160] 

method 

heat  capacity  are  based  on 
previous  time  step  (i.e., 
explicit  scheme) 

Modified 

FDM:  ID 

Heat  capacity 

Implicit 

•  Developed  for  Internal 

Experimental 

[161] 

“TYPE101” 

method 

partition  wall 

“TYPE1270” 

Lumped 

Quasi-heat  source 

•  Very  simplified  model 

N/A 

[162] 

method 

method 

•  Internal  layer  within  an 

using  heat 

envelope 

balance 

•  Based  on  lumped  heat 

balance  (not  a  finite  volume), 
low  accuracy 

•  For  phase  change  at  fixed 

temperature 

ESP-r 

SPMCMP53- 

FDM:  ID 

Heat  capacity  and 

•  Low  time  step 

N/A 

[164,165] 

SPMCMP56 

heat  source 
method 

BSim 

FVM:  ID 

Heat  capacity 

Implicit 

•  Low  time  step  to  avoid 

Experimental 

[171] 

method 

instability 

RADCOOL 

FDM:  ID 

Heat  capacity 
method 

Implicit 

[172] 

ESim 

FDM:  ID 

Heat  capacity 

Explicit 

•  Explicit  scheme  requires  low 

Experimental 

[179] 

method 

time  step  to  avoid  instability 

FVM:  finite  volume  method,  FDM:  finite  difference  method. 
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literature  available  for  these  programs,  only  few  will  be  discussed 
hereafter. 


5.5.1.  RADCOOL 

RADCOOL  is  a  design  tool  for  cooling  and  heating  system 
developed  at  the  Lawrence  Berkeley  National  Laboratory  of  the  US 
[175].  The  program  was  created  using  the  simulation  problem 
analysis  and  research  kernel  (SPARK)  [176].  A  one  dimensional 
finite-difference  model  for  a  wall  with  PCM  was  added  to  this 
validated  thermal  building  simulation  program  [172].  The  model 
was  then  used  to  study  the  capability  of  a  double  PCM-wallboard 
to  achieve  thermal  comfort  without  using  mechanical  cooling 
system  under  a  typical  climatic  condition  of  Sunnyvale,  California 
[177], 

5.5.2.  ESim 

ESim  was  developed  at  University  of  Dayton  for  building 
energy  simulation  and  can  be  downloaded  from  the  developer 
website  [178].  The  simulation  program  was  expanded  and  vali¬ 
dated  to  model  PCM-wallboards  using  an  explicit  finite-difference 
approach  [179].  A  list  of  template  files  are  available  for  use  but 
with  limited  capability  to  model  complex  buildings  and  systems. 

5.5.3.  PCMExpress 

PCMExpress  is  a  planning  and  simulation  program  for  build¬ 
ings  using  phase  change  materials  [174],  which  was  developed  by 
a  German  company,  in  collaboration  with  the  Fraunhofer  Institute 
for  Solar  Energy  (ISE)  in  Freiburg  and  partners  from  industry 
[180],  The  program  simulates  a  free  floating  building  with  a 
library  data  for  weather  and  various  construction  materials 
including  PCM  and  the  flexibility  to  add  new  materials.  It  is  an 
effective  tool  to  evaluate  the  economic  and  technical  feasibility  of 
PCM  usage  during  an  early  design  stage.  The  mathematical  model 
of  the  heat  transfer  process  in  PCM  is  not  available.  The  model  has 
been  tested  by  Castell  [141]  who  found  that  the  simulations 
deviate  significantly  from  the  experiments.  As  commented  by 
Castell,  the  discrepancy  could  be  attributed  to  the  lack  of  accurate 
infiltration  model  in  the  program.  The  program  however  has  been 
used  to  demonstrate  the  impact  of  using  PCM  in  residential  and 
commercial  buildings  [180,181]. 

5.6.  Summary 

Whole  building  simulation  programs  play  an  important  role 
for  studying  the  economic  and  technical  feasibility  of  PCMs.  This 
section  reviews  the  capability  of  various  whole  building  simula¬ 
tion  programs  as  summarized  in  Table  3.  It  is  noted  that  most 
PCM  models  integrated  into  whole  building  simulation  programs 
are  based  on  the  heat  capacity  method.  Hence,  it  becomes 
necessary  to  reduce  the  typical  one  hour  time  step  to  a  very 
small  time  step  (i.e.,  in  order  of  minutes)  to  achieve  acceptable 
level  of  accuracy.  For  one  year  thermal  performance  evaluation, 
building  simulation  programs  become  computationally  inefficient 
since  iterative  methods  are  used  in  each  time  step.  Additionally, 
the  convergence  may  not  be  achieved  due  to  numerical  instability 
especially  when  PCM  enters  or  leaves  the  phase  change  region. 
With  all  constraints  and  limitations  above,  none  of  whole  building 
simulation  programs  are  currently  implementing  efficient  math¬ 
ematical  models  that  are  quick,  accurate  and  numerically  stable 
at  realistic  time  step.  It  becomes  important  to  thoroughly  inves¬ 
tigate  different  mathematical  models  with  various  numerical 
approaches  for  modeling  PCMs  in  whole  building  simulation 
programs. 


6.  Conclusions 

Significant  heat  storage  offered  by  phase  change  materials  is 
promising  and  favorable  for  developing  various  innovative  building 
energy  technologies.  To  quantify  the  technical  and  economic  feasi¬ 
bility  of  PCM-embedded  technologies,  it  requires  the  development 
of  proper  computational  models.  This  study  reviews  various  numer¬ 
ical  modeling  approaches  of  phase  change  problems  such  as  the 
enthalpy,  heat  capacity,  temperature-transforming  and  heat  source 
methods.  The  main  features,  advantages  and  disadvantages  of  each 
method  have  been  discussed.  The  discretized  form  of  the  heat 
equation  with  PCM  can  either  be  solved  with  nonlinear  solvers  such 
as  Newton’s  methods  or  via  linearizing  the  nonlinear  term  and  using 
linear  solvers  such  as  iterative  methods.  For  both  approaches,  the 
numerical  solutions  are  computationally  inefficient  or  difficult  to 
reach  convergence.  Therefore,  fast  numerical  schemes  are  suggested 
such  as  the  quasi  enthalpy  non-iterative  scheme  or  the  enthalpy 
conservative  iterative  scheme. 

Using  these  general  mathematical  methods,  different  compu¬ 
tational  models  have  been  developed  to  simulate  PCMs  in  build¬ 
ing  enclosures.  Based  on  the  level  of  complexity,  models  are 
classified  into  three  categories:  simple,  intermediate  and  sophis¬ 
ticated  models.  Majority  of  these  models  have  been  validated 
using  analytical  solutions,  comparative  testing  using  validated 
numerical  models,  and/or  experimental  results. 

While  many  models  are  used  to  study  the  heat  transfer  in  an 
enclosure  unit,  a  few  models  have  been  integrated  into  whole 
building  simulation  programs.  A  variety  of  models  are  available 
and  some  are  available  with  no  cost  to  users  such  as  EnergyPlus, 
“TYPE204”  in  TRNSYS  or  ESP-r.  These  particular  models,  however, 
have  limitations  on  modeling  PCM  including  the  time  and  spatial 
resolutions,  inability  to  model  hysteresis,  lack  of  validations  of 
some  models,  and  poor  computational  efficiency.  These  modeling 
challenges  add  complexity  to  the  already  existing  uncertainties  in 
experimental  results  of  PCM’s  thermal  behaviors.  Therefore, 
further  research  is  needed  to  quantitatively  explore  the  prediction 
performance  of  different  models  including  their  limitations  on 
accuracy,  parameters  sensitivity,  speed,  and  stability  for  modeling 
PCM  envelopes  under  different  climatic  and  operating  conditions. 
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